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Acoustic  Observation  of  the  Time  Dependence 
of  the  Roughness  of  Sandy  Seafloors 

Darrell  R.  Jackson,  Michael  D.  Richardson,  Kevin  L.  Williams,  Anthony  R  Lyons,  Member,  IEEE, 
Christopher  D.  Jones,  Kevin  B.  Briggs,  and  Dajun  Tang 


Abstract — A  statistical  model  for  the  time  evolution  of  seafloor 
roughness  due  to  biological  activity  is  applied  to  photographic 
and  acoustic  data.  In  this  model,  the  function  describing  small 
scale  seafloor  topography  obeys  a  time-evolution  equation  with 
a  random  forcing  term  that  creates  roughness  and  a  diffusion 
term  that  degrades  roughness.  When  compared  to  acoustic  data 
from  the  1999  and  2004  Sediment  Acoustics  Experiments  (SAX99 
and  SAX04),  the  model  yields  diffusivitics  in  the  range  from 
3.5x10“^^  to  2.5  X  10~^°  m^s“^  (from  10  to  80  cm^yr“D»with 
the  larger  values  occurring  at  sites  where  bottom-feeding  fish  were 
active.  While  the  experimental  results  lend  support  to  the  model, 
a  more  focused  experimental  and  simulation  effort  is  required  to 
test  several  assumptions  intrinsic  to  the  model. 

Index  Terms — Acoustic  scattering,  bioturbation,  diffusion  pro¬ 
cesses,  seafloor,  sediment  transport. 


I.  Introduction 

Biological  activity  (bioturbation)  at  the  seafloor 
causes  small  scale,  random  sediment  transport  that  is 
often  approximated  as  a  diffusive  process.  In  the  absence 
of  hydrodynamic  forcing  by  storms  or  strong  tidal  currents, 
biological  processes  will  dominate  the  evolution  of  bottom 
roughness  features.  As  noted  in  [1],  this  was  the  case  at  one  of 
the  sites  considered  in  this  paper.  The  evolution  of  roughness 
has  been  treated  as  a  horizontal  diffusive  process  [2],  [3].  In 
this  view,  horizontal  diffusion  of  sediment  particles  due  to 
the  activity  of  small  animals  living  within  or  near  the  seafloor 
will  tend  to  reduce  roughness  features.  Counteracting  this 
degradation  of  roughness,  biological  processes  also  create  new 
roughness  features.  If  both  processes  act  for  a  sufficiently  long 
time  and  at  constant  rates,  a  statistical  equilibrium  will  be 
reached  in  which  details  of  roughness  change  with  time  while 
the  power  spectral  density  of  random  roughness  is  constant  [1]. 
Thi.s  assumption  is  used  in  [4]  where  relations  between  acoustic 
and  roughness  temporal  correlations  are  developed  and  applied. 
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Fig.  1.  Pockmarks  created  at  the  sandy  SAX99  site  by  feeding  pinfish.  Note 
the  remnant  of  quasi-periodic,  ripple-like  structure  that  was  obliterated  within 
24  hours  by  this  biological  activity  [6]. 


In  this  paper,  similar  relations  are  used  in  conjunction  with  a 
model  for  the  time  evolution  of  seafloor  roughness  [5].  This 
model  will  be  examined  in  some  detail  and  compared  with 
photographic  and  acoustic  data. 

A  photographic  example  of  biological  activity  creating 
and  modifying  seafloor  roughness  is  shown  in  Fig.  1.  In  the 
1999  Sediment  Acoustics  Experiment  (SAX99),  conducted  in 
shallow  water  in  the  northeastern  Gulf  of  Mexico,  fish  feeding 
on  small  benthic  organisms  created  roughness  in  the  form  of 
overlapping  pockmarks. 

Numerous  studies  have  interpreted  and  quantified  vertical 
mixing  in  sediments  as  a  diffusive  process,  but  horizontal 
diffusion  is  more  relevant  to  the  problem  of  roughness  time 
evolution.  Wheatcroft  [3]  has  presented  data  suggesting  that 
horizontal  diffusivities  are  larger  than  vertical  diffusivities  by 
an  approximate  factor  of  ten.  While  most  studies  of  sediment 
diffusion  have  focused  on  fine-grained  clays  and  muds,  this 
paper  is  concerned  with  sandy  sediments.  Of  direct  relevance 
to  this  work  are  measurements  of  horizontal  diffusivity  at 
the  2004  Sediment  Acoustics  Experiment  (SAX04)  site  [2]. 
These  measurements  used  high-frequency  sonar  images  of 
time-varying  roughness  and  obtained  diffusivities  tn  the  range 
from  1  to  20x10^^^^  m^s”^  (from  300  to  6000  cm^yr“^). 
The  same  data  set  is  discussed  extensively  in  [1]  although 
diffusivity  values  are  not  given.  Applying  a  method  developed 
below,  a  diffusivity  of  3.2  x  10”‘‘^  m^s“^  (1000  cm^yr“*)  can 
be  obtained  from  the  spectral  decay  data  of  [1].  In  this  paper. 
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diffusivity  will  be  denoted  D  and  given  in  SI  units.  The  more 
traditional  units  (cm^/yr)  will  also  be  given  occasionally  to 
facilitate  comparison  with  the  literature. 

Mathematical  models  have  been  developed  to  describe  trans¬ 
port  of  sediment  particles  by  bioturbation.  The  most  common 
are  1-D  biodiffusion  transport  models  that  employ  a  depth-de- 
pendent  biodiffusion  or  mixing  coefficient  which  characterizes 
the  intensity  of  vertical  mixing  as  a  function  of  depth  [8]. 
Most  models  for  biodiffusion  are  based  on  measurements  of 
particle  transport  in  the  vertical  direction,  using  solid-phase 
continuous  or  impulsive  tracers  such  as  naturally  occurring 
radionuclides  (Pb-210,  Th-235,  Cs-137,  Be-7),  radionuclides 
associated  with  fallout  from  nuclear  testing  (Pu>239,  Pu’24()), 
or  by  introduced  tracers  such  as  luminophores  or  labeled 
glass  beads.  Vertical  mixing  coefficients  measured  in  coastal 
sediments  typically  range  between  3  and  30xl0~^^  m^s~^ 
(from  10  to  100  cm^yr~^)  over  characteristic  depth  scales  of 
10-30  cm  [9]-[ll].  This  rapid  mixing  of  surficial  sediments 
tends  to  obliterate  signatures  of  depositional  layers  over  scales 
of  tens  of  years.  Measured  diffusivities  often  decrease  with 
depth  in  the  sediment  and  can  be  modeled  as  decreasing  or 
layer-specific  mixing  coefficients. 

A  seafloor-mounted,  rotating  sonar  platform  operating  at 
40  kHz  was  used  to  observe  changes  in  seafloor  backscatter 
at  the  Sediment  Transport  on  Shelves  and  Slopes  (STRESS) 
site  off  the  northern  California  coast  [12],  [13].  The  complex 
correlation  between  seafloor  echoes  acquired  at  different  times 
(defined  later  in  this  paper  and  referred  to  as  the  “ping-to-ping 
correlation”)  yielded  a  decorrelation  time  on  the  order  of  one 
month.  Decorrelation  times  of  the  same  order  were  reported 
for  a  site  off  Orcas  Island  in  Washington  State  [5],  [14],  [15]. 
At  both  the  STRESS  and  Orcas  sites,  the  seafloor  was  a  soft 
mud,  and  the  changes  were  attributed  to  biological  activity 
altering  the  details  of  volume  scattering  within  the  sediment. 
As  the  present  interest  is  in  changes  in  seafloor  roughness, 
measurements  on  harder  sandy  seafloors  are  of  primary  in¬ 
terest.  For  such  seafloors,  acoustic  scattering  due  to  roughness 
is  expected  to  dominate  volume  scattering  for  grazing  angles 
smaller  than  the  critical  angle  (about  30°  for  the  sandy  sites 
considered  here).  Short-term  backscatter  fluctuations  at  two 
sandy  sites  were  ascribed  by  Stanic  and  Kennedy  [16],  [17]  to 
water-column  fluctuations  rather  than  changes  in  the  seafloor. 
As  part  of  the  coastal  benthic  boundary  layer  (CBBL)  pro¬ 
gram,  longer  term  observations  at  a  sandy  site  were  carried 
out  at  40  kHz  [18].  The  results  were  presented  in  terms  of 
ping-to-ping  correlation  and  exhibited  a  decorrelation  time 
on  the  order  of  one  day.  Numerous  shell  pieces  lay  upon  the 
seafloor  and  it  is  possible  that  these  discrete  scatterers  played 
a  significant  role  in  scattering.  During  SAX99,  similar  corre¬ 
lation  measurements  were  made  and  are  reported  here.  Also 
during  SAX99,  the  seafloor  was  roughened  artificially  and  the 
resulting  time-dependent  backscattering  strength  was  reported 
[6].  These  data  and  previously  unpublished  data  from  SAX04 
will  be  compared  with  the  model  of  Jones  [5]. 

Photographic  time  series  are  also  relevant  to  the  present  sub¬ 
ject.  A  series  of  digital  stereophotographs  were  obtained  by 
Lyons  during  SAX99  and  the  corresponding  time-dependent 
roughness  spectra  were  reported  in  [6].  A  longer  time  series 
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of  digital  stereophotgraphs  has  been  obtained  at  a  sandy  site 
near  the  island  of  Elba,  Italy  [19],  and  at  the  SAX04  site  [4].  In 
[19],  changes  in  roughness  spectra  were  attributed  to  hydrody¬ 
namic  activity  (forming  ripples)  during  storms  and  bioturbation 
(mound  building  by  mud  shrimp)  during  quiet  periods.  It  was 
suggested  that  the  sediment  mixing  is  primarily  vertical  rather 
than  horizontal,  as  the  mud  shrimp  create  funnel-shaped  features 
and  mounds  from  U-shaped  burrows  while  feeding.  Temporal 
changes  in  backscatter  strength  were  calculated  using  a  model 
with  the  time-varying  measured  spectra  as  input.  The  changes 
in  root  mean  square  (RMS)  roughness  were  found  to  be  small 
compared  to  diurnal  changes  in  parameters  describing  the  shape 
of  the  power  spectral  density  function. 

Section  II  defines  the  model  in  two  parts:  the  first  part  is 
the  roughness  evolution  model  of  [5],  and  the  second  part 
gives  the  connection  between  roughness  and  acoustic  backscat¬ 
tering  via  conventional  small-roughness  perturbation  theory. 
Section  II-A  considers  the  implications  of  the  model  for  the 
evolution  of  artificially  created  roughness,  Section  II-B  does 
the  same  for  naturally  created  roughness,  and  Section  II-C 
examines  the  statistics  of  the  forcing  term  in  the  equation  for 
time  evolution.  Section  III  treats  issues  related  to  acoustic 
measurement  of  time-dependent  backscattering,  and  Section  IV 
gives  brief  descriptions  of  the  acoustic  measurement  systems. 
Properties  of  the  experimental  sites  relevant  to  model-data 
comparisons  are  discussed  in  Section  V,  and  these  comparisons 
are  made  in  Section  VI  for  artificially  roughened  seafloors  and 
in  Section  VII  for  naturally  evolving  roughness.  Discussion  of 
the  present  status  of  the  model  and  suggestions  for  future  work 
are  given  in  Section  VIII. 

II.  Model  Definition 

The  model  of  Jones  [5]  allows  for  a  position-dependent  dif¬ 
fusivity,  and  the  forcing  term  is  written  as  a  superposition  of  el¬ 
emental  “shape  functions”  incorporating  the  microtopography 
of  individual  roughness-creating  events.  A  simplified  version 
will  be  defined  here,  with  position-independent  diffusivity  and 
with  the  forcing  term  defined  in  rather  general  terms,  but  con¬ 
strained  by  assumptions  regarding  its  independence  relative  to 
the  existing  relief.  Time-dependent  seafloor  roughness  will  be 
described  by  the  function  /(R,  t)  giving  the  vertical  displace¬ 
ment  of  the  seafloor  from  the  mean  horizontal  plane.  The  vari¬ 
able  R  is  a  2-D  vector  composed  of  the  x-  and  y-coordi nates 
(horizontal),  R  =  (x*,  y).  Throughout  the  remainder  of  this 
paper,  /(R,  t)  will  be  referred  to  as  the  “relief  function.” 

The  model  is  defined  by  the  differential  equation  governing 
the  time  dependence  of  the  relief  function  as  well  as  by  the  as¬ 
sumed  statistical  properties  of  the  terms  in  this  equation.  Spe¬ 
cializing  an  equation  given  in  [5],  the  2-D  relief  function  will  be 
assumed  to  obey 

a/(R.t)  ^^vV(R,0  +  5(R./.).  (1) 

The  first  term  on  the  right-hand  side  in  (I),  the  diffusion  term, 
causes  roughness  to  decay  with  time,  while  the  second  term,  the 
forcing  term,  counteracts  this  decay.  The  diffusivity  D  is  the 
defining  parameter  of  the  model  and  is  also  the  parameter  to 
be  estimated  in  fits  of  the  model  to  data.  As  both  diffusion  and 
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forcing  are  assumed  to  be  due  to  biological  activity,  it  is  rea¬ 
sonable  to  inquire  as  to  which  sort  of  activity  leads  to  diffusion 
and  which  leads  to  roughness  generation.  No  definitive  answer 
is  offered  here,  but  Hay  [1]  ascribes  ripple  decay  in  his  SAX04 
data  to  bottom-feeding  fish:  In  [14],  it  is  suggested  that  diffusion 
results  from  the  activity  of  smaller  animals  and  mechanical  col¬ 
lapse  of  features  created  by  larger  animals,  while  forcing  rep¬ 
resents  the  activity  of  larger  animals.  In  the  present  problem, 
it  can  be  asserted  that  the  forcing  term  represents  the  creation 
of  feeding  pits  and  other  features.  The  forcing  term  will  be  as¬ 
sumed  to  be  statistically  uncorrelated  with  the  relief  function 

(/(Ri,05(R2,  0>  =  0.  (2) 

The  symbols  ()  denote  a  formal  average  of  an  ensemble  of  dif¬ 
ferent  random  realizations  of  the  relief  function  and  forcing 
term.  It  is  important  that  both  are  evaluated  at  equal  times  t  in 
(2).  There  will  be  some  correlation  between  the  forcing  term  at 
any  given  time  and  the  relief  at  a  later  time,  as  part  of  this  relief 
results  from  the  forcing.  Equations  (1)  and  (2)  embody  assump¬ 
tions  whose  primary  motivation  is  simplicity,  but  which  may 
prove  to  be  approximately  true  in  real  situations.  First,  (1)  says 
that  the  effect  of  the  forcing  term  is  simply  additive.  In  combi¬ 
nation,  (1)  and  (2)  require  that  feeding  pits,  etc.,  cannot  be  pref¬ 
erentially  located  with  respect  to  features  of  the  topography.  For 
example,  it  is  conceivable  that  feeding  by  fish  as  in  Fig.  1  will 
take  place  preferentially  at  highs  or  lows  of  the  topography,  but 
this  is  not  allowed  by  the  model.  Equations  (1)  and  (2)  also  do 
not  allow  the  shape  of  any  feeding  feature  (such  as  a  pit)  to  be 
influenced  by  the  preexisting  topography.  That  is,  the  shape  of 
the  pit  (defined  by  the  change  in  2:-coordinate  over  the  pit  with 
respect  to  the  preexisting  topography)  is  the  same  whether  the 
pit  occurs  at  a  topographic  high,  a  topographic  low,  or  on  a  slope. 
While  it  is  certain  that  these  assumptions  cannot  be  strictly  true, 
the  issue  is  whether  they  provide  useful  approximations. 

The  final  assumption  is  that  the  forcing  term  5(R,  t)  is  a 
stationary  random  process  with  respect  to  spatial  coordinates 
and  time.  That  is,  biological  creation  of  roughness  proceeds  at 
a  steady  pace  for  all  times  of  interest  and  has  constant  strength 
over  the  seafloor  area  of  interest.  While  the  assumption  of  tem¬ 
poral  stationarity  is  also  used  in  [4],  those  authors  note  that  it 
was  not  well  satisfied  in  their  SAX04  data.  Temporal  changes 
in  roughness  spectra  have  also  been  reported  in  [19]  and  [20]. 
On  the  other  hand.  Hay  [1]  finds  that  roughness  spectra  mea¬ 
sured  during  SAX04  exhibit  stationarity  in  the  high-frequency 
band  (8-24  cycles/m,  this  is  wave  number  divided  by  27r).  In  this 
paper,  we  depart  from  SI  practice  in  denoting  the  units  of  spa¬ 
tial  frequency  by  either  cycles  per  meter  or  radians  per  meter,  as 
appropriate.  A  major  cause  of  temporal  nonstationarity  is  the  di¬ 
urnal  behavior  of  animals  living  near  or  within  the  seafloor.  For 
example,  the  diurnal  changes  in  roughness  documented  in  [19] 
were  due  to  changes  in  the  activity  of  burrowing  shrimp.  Fish 
commonly  exhibit  diurnal  behavior,  forming  schools  during  the 
day  and  descending  in  the  water  column.  Also,  small  animals 
are  known  to  emerge  from  the  seafloor  into  the  water  column 
at  night  [21].  If  roughness  statistics  are  observed  and  averaged 
over  several  days,  it  may  be  reasonable  to  treat  animal  activity 
as  stationary  on  these  time  scales.  In  this  paper,  the  stationarity 


assumption  is  made  primarily  on  the  grounds  of  simplicity  and 
could  be  avoided  by  allowing  the  forcing  term  to  have  temporal 
changes  in  its  statistics  on  time  scales  much  longer  than  those 
characterizing  the  roughness-causing  events.  As  these  events 
have  time  scales  measured  in  seconds  or  less,  this  approach  is 
feasible,  but  not  warranted  at  present  owing  to  lack  of  informa¬ 
tion  on  the  slow  time  dependence  of  biological  activity  at  the 
experimental  sites. 

When  artificial  roughness  is  introduced,  it  will  be  modeled 
as  an  additive  initial  term  in  the  relief  function  and  not  as  a 
term  in  the  forcing  function.  If  there  is  no  artificial  modification 
of  roughness,  and  if  the  forcing  term  is  statistically  stationary 
and  has  acted  for  a  sufficiently  long  time,  the  spectrum  of  nat¬ 
urally  created  roughness  will  reach  a  time-independent  equilib¬ 
rium  value.  The  relief  function  itself  will  change  with  time,  but 
the  resulting  spectrum  will  be  constant  in  time.  From  the  point 
of  view  of  the  present  development,  this  equilibrium  spectrum 
can  be  considered  a  parameter  of  the  model  because  it  will  be 
used  to  determine  the  statistics  of  the  forcing  term  (actually,  the 
spectrum  provides  many  parameters,  as  it  is  always  measured 
over  a  range  of  spatial  frequencies). 

The  link  to  acoustic  scattering  is  provided  by  small-roughness 
perturbation  theory,  which  can  be  used  to  show  that  the  complex 
backscattered  pressure  at  time  t  is  of  the  form 

(3) 

In  this  expression,  //  is  a  factor  depending  on  the  acoustic  fre¬ 
quency,  grazing  angle,  and  acoustic  parameters  of  the  water  and 
sediment  (mass  density,  sound  speed,  and  attenuation).  A  nondi- 
rectional,  continuous- wave  source  is  assumed  with  slant  range  ? 
from  the  source  to  the  scattering  location  of  interest.  The  Fourier 
transform  of  the  relief  function  is  represented  by  the  notation 
F(K,  i),  where 

F{K,  t)  =  I  /(R,  (4) 

and  where  K  =  {Kx^  I^y)  is  a  2-D  wave  vector.  In  (3),  the 
Fourier  transform  is  evaluated  at  the  Bragg  wave  vector  K/j  = 
(2A:cos^,0),  where  k  =  uj/c  is  the  acoustic  wave  number  in 
water,  and  0  is  the  grazing  angle.  Equation  (3)  is  derived  under 
the  assumption  that  the  relief  function  falls  to  zero  outside  a 
small  seafloor  patch  situated  at  slant  range  r  from  the  source. 
This  contradicts  the  earlier  assumption  that  the  relief  function 
is  a  stationary  random  process  with  respect  to  the  horizontal 
spatial  coordinates.  This  contradiction  can  be  avoided  by  as¬ 
suming  a  directional  source,  but  the  formalism  becomes  more 
involved.  In  the  present  approach,  one  can  take  the  relief  func¬ 
tion  to  be  the  product  of  a  stationary  process  and  a  smooth  spa¬ 
tial  window  function.  If  the  window  is  wide  compared  to  the 
horizontal  scales  of  roughness  of  interest,  yet  small  compared 
to  the  slant  range,  the  assumption  of  spatial  stationarity  can  still 
be  employed. 

From  (3),  the  covariance  of  backscattered  pressure  {p2P\) 
obtained  at  times  t\  and  t2  is  of  the  form 

{P2Pi)  =  QW{KB..tut2)  (5) 
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where  Q  is  an  inessential  factor  depending  upon  grazing  angle, 
frequency,  and  acoustic  properties.  The  cross  spectrum  of  relief 
in  (5)  is  related  to  the  Fourier  transform  of  the  relief  function  as 
follows: 


{F(K2,  t2)F*{Ku  h))  =  W{Ku  tu  t2)S{Ki-K2)  .  (6) 

The  presence  of  the  Dirac  delta  function  is  only  appropriate  if 
/(R,  is  stationary  with  infinite  extent  in  the  spatial  coordi¬ 
nates.  As  a  result  of  the  windowing  assumed  here,  the  delta  func¬ 
tion  is  actually  an  approximation  to  a  sharply  peaked,  but  finite, 
function.  The  cross  spectrum  is  related  to  the  covariance  of  the 
relief  function  as  follows: 

ir(K.  ti,  t2)  =  /b(R,  (7) 

(27r )  J 

where  K  =  /v  ,J  is  a  2-D  wave  vector,  and  i?(R,  ^1,^2) 
is  the  two-time  covariance  of  the  relief  function 

/J(R2  -  Ri-  /-i,  ^2)  =  (/(R2,  ^2)/(Ri,  ^i))  .  (8) 

As  a  result  of  the  stationarity  approximation,  only  the  coordi¬ 
nate  difference  appears  on  the  left-hand  side  of  (8).  The  time 
dependence  is  left  general  to  accommodate  cases  in  which  arti¬ 
ficial  roughness  is  introduced,  causing  a  breakdown  of  temporal 
stationarity. 

The  expressions  given  above  will  be  used  to  treat  two  dif¬ 
ferent  problems  in  time  dependence  of  backscattering.  In  the 
first  problem,  artificial  roughness  is  deliberately  created,  and 
the  time  decay  of  this  roughness  is  observed  as  a  decay  in  the 
backscattering  cross  section  to  be  defined  later.  The  scat¬ 
tering  cross  section  is  proportional  to  the  equal-time  spectrum 
obtained  by  setting  =  ^2  in  (7).  This  spectrum  is  evalu¬ 
ated  at  the  Bragg  wave  vector  defined  earlier.  In  the  second 
problem,  the  time  dependence  of  natural  roughness  is  observed 
through  change  in  the  details  of  the  seafloor  echo  in  successive 
pings.  These  changes  are  quantified  by  means  of  the  complex 
ping-to-ping  correlation  coefficient,  proportional  to  the  cross 
spectrum  (7)  with  ti  7^  ^2-  B  will  be  possible  to  discuss  both 
of  these  problems,  even  comparing  the  model  with  measure¬ 
ment  data,  without  examining  the  forcing  term.  Later,  it  will  be 
shown  that  the  assumptions  embodied  in  (1)  and  (2)  suffice  to 
determine  the  statistical  properties  of  the  forcing  term  (provided 
the  equilibrium  spectrum  is  known).  Using  this  idea,  the  forcing 
term  will  be  discussed  in  some  detail. 

A.  Time  Decay  of  Artificial  Roughness 

During  SAX99  [22],  [23]  and  SAX04  [24],  [25],  artificial 
roughness  was  introduced  by  raking  the  seafloor  to  create  reg¬ 
ular  ripple  patterns  that  produced  strong  backscattering  at  the 
frequency  corresponding  to  the  Bragg  wavelength  [6],  [25].  To 
address  this  problem  without  direct  appeal  to  the  forcing  term, 
the  relief  function  will  be  expressed  as  the  sum  of  natural  and 
artificial  terms 


/(R,  0  =  /n(R.  f)  +  /a(R,  0.  (9) 


The  artificial  roughness /„(R,  is  created  at  ^  =  0,  and  the  nat¬ 
ural  roughness  /n(R,  0  is  assumed  to  evolve  (even  through  the 
moment  of  artificial  roughness  creation)  as  it  would  normally. 
The  superposition  expressed  in  (9)  is  questionable  for  small 
times  because  the  creation  of  artificial  roughness  will  modify 
the  existing  natural  roughness.  This  idealization  is  reasonable, 
however,  if  the  artificial  roughness  is  initially  greater  than  the 
natural  roughness.  In  this  case,  the  time  required  for  the  artificial 
roughness  to  decay  to  a  value  comparable  to  the  natural  rough¬ 
ness  is  greater  than  the  time  required  for  the  natural  roughness 
to  regain  statistical  equilibrium.  In  other  words,  by  the  time  the 
artificial  roughness  decays  to  a  level  where  it  no  longer  masks 
the  natural  roughness,  the  natural  roughness  will  have  returned 
to  equilibrium.  The  natural  roughness  is  time  dependent  owing 
to  the  competing  processes  of  generation  and  diffusion,  but  it 
has  a  time-independent  power  spectral  density  II^ti(K).  As  the 
natural  roughness  term  in  (9)  satisfies  (1),  the  artificial  rough¬ 
ness  must  obey  the  diffusion  equation  with  no  forcing 

=  (10) 

It  will  be  assumed  that  the  two  roughness  components  are  un¬ 
correlated 


(/n(Ro  +  Ri  fl)/a(Ro?  ^2))  —  fl  (11) 

insuring  that  the  overall  roughness  spectrum  is  the  sum  of  the 
spectra  for  /n(R,  t)  and  /„(R,  t) 

IV(K,  L  t)  =  iy„(K)  +  VT„(K,  ft),  (12) 

Using  (6),  the  time-dependent  spectrum  L  0  can  be 

written  in  terms  of  the  Fourier  transform  of  /a(R,  t) 


(F„(K.  0^«*(Ko,  t))  =  VF«(K,  L  O^(K-Ko).  (13) 

The  diffusion  equation  (10)  gives 

Fa(K,  t)  =  F„(K,  0)exp(-/^2^<)-  (14) 

The  final  result  for  the  time  dependence  of  the  roughness  spec¬ 
trum  is 


H/(K,  t,  t)  =  H'„(K)  +  Wa{K.  0,  0)cxp{-2K^Dt)  (15) 


where  iy„.(K,  0,  0)  is  the  spectrum  of  the  artificial  roughness 
immediately  after  its  creation.  Note  that  the  1/e  decay  time  is 


1 

~  2K^D' 


(16) 


In  practice,  the  time-dependent  backscattered  acoustic  intensity 
is  measured,  but  it  is  simply  proportional  to  W(fK,  f,  f)  evalu¬ 
ated  at  the  Bragg  wave  vector  K^.  Specializing  (5)  to  the  present 
case 


{PP^)  =  QW{K^,ft).  (17) 

Taken  together,  (15)  and  (17)  show  that  the  mean  square 
backscattered  pressure  is  the  sum  of  a  constant  term  due  to 
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natural  roughness  and  a  decaying  term  due  to  artificial  rough¬ 
ness.  The  initial  value  of  the  decaying  term  is  determined  by 
the  initial  artificial  roughness.  The  backscattering  cross  section 
(which  is  proportional  to  {PP*))  is  dependent  on  the  grazing 
angle  0,  and  is  defined  through  the  relation  [7] 


a{e)  = 


{PP')r^ 

~mF 


(18) 


where  Pt  is  the  complex  amplitude  of  the  pressure  incident  upon 
a  scattering  patch  of  area  A  situated  at  range  r  from  the  receiver 
at  which  P  is  measured.  In  the  case  at  hand,  the  cross  section 
has  the  form 


<r(0)  =  a„{0)+<Ta{e)e.xp{-2K^Dt)  (19) 

with  a  constant  “natural”  term  and  a  decaying  “artificial” 
term.  Equation  (15)  predicts  that  short-wavelength  (large  K) 
roughness  features  will  decay  more  rapidly  than  long- wave¬ 
length  (small  K)  features.  Since  the  Bragg  wave  number  is 
proportional  to  acoustic  frequency,  this  leads  to  a  prediction  for 
the  frequency  dependence  of  the  decay  rate  of  scattering  due 
to  artificially  created  roughness.  According  to  (19),  the  decay 
rate  should  be  proportional  to  the  square  of  acoustic  frequency. 
The  application  of  (19)  to  backscatter  data  will  be  discussed  in 
Section  VI. 


B.  Time  Evolution  of  Natural  Roughness 

In  the  natural  state,  the  seafloor  is  not  disturbed  by  the  ex¬ 
perimentalist,  and  changes  in  the  details  of  seafloor  roughness 
can  be  observed  by  means  of  correlation  between  sonar  echoes 
taken  at  different  times.  This  is  the  ping-to-ping  correlation,  for 
which  the  complex  cross  correlation  of  backscattered  pressure 
is  of  interest.  In  small-roughness  perturbation  theory,  this  cor¬ 
relation  can  be  expressed  in  terms  of  the  two-time  cross  spec¬ 
trum  using  (5).  Owing  to  the  assumed  temporal  stationarity  of 
natural  roughness,  the  cross  spectrum  will  only  depend  on  the 
“lag,”r  =  ^2  “  f  1  •  The  problem  of  measuring  the  pressure  cross 
correlation  will  be  discussed  in  Section  III.  The  essential  factor 
is  the  lag-dependent  cross  spectrum,  which  can  be  obtained  from 
the  diffusion  equation  by  separating  the  relief  function  into  a  de¬ 
caying  part  that  is  due  to  all  forcing  occurring  for  negative  times 
and  a  growing  part  that  is  due  to  all  subsequent  forcing 

/(R,  t)  =  /d(R,  t)  +  /^(R,  t).  (20) 

The  assumed  lack  of  correlation  between  the  forcing  term  and 
the  relief  function  at  any  given  time  insures  that  the  two  compo¬ 
nents  of  the  relief  function  are  uncorrelated.  As  a  result  of  this 
lack  of  correlation,  and  because  /y(R,  0)  =  0,  the  cross  spec¬ 
trum  for  times  0  and  t  can  be  obtained  from  (6)  using 


(F(K2.  0))  =  (Frf(K2.  <)F;(K,,  0)) .  (21) 

By  definition,  the  decaying  part  of  the  roughness  obeys  the  dif¬ 
fusion  equation  with  no  forcing  term  for  i  >  0,  so  F^i(K,  t)  is 
of  the  form  (14).  The  result  for  the  cross  spectrum  is 


thus 


(F2A*)  =  -  (1)] .  (23) 

Analogous  to  the  previous  result  for  the  decay  of  scattering  due 
to  artificially  created  roughness,  (23)  gives  a  prediction  for  the 
frequency  dependence  of  ping-to-ping  correlation.  Considering 
the  two  frequencies  used  in  the  experiments  to  be  reported  later, 
ping-to-ping  correlation  should  decay  56  times  faster  at  300  kHz 
than  at  40  kHz. 


C.  The  Forcing  Term 

The  forcing  term  5(R,  t)  in  (1)  represents  new  roughness 
features  created  near  time  t.  Specifically,  2:  =  5(R,  t)Af  is  the 
relief  of  the  mounds,  pits,  surface-penetrating  tubes,  etc.,  cre¬ 
ated  in  the  time  interval  between  t  and  1 4-  At.  A  detailed  model 
for  the  time  dependence  of  roughness  might  specify  the  forcing 
term  as  a  superposition  of  elemental  shapes  representing  these 
features  as  in  [5].  In  this  work,  such  details  will  not  be  consid¬ 
ered,  but  the  spectrum  of  the  forcing  term  can  be  determined 
and  may  offer  clues  as  to  the  nature  of  the  forcing. 

Consider  the  growing  part  /<;(R,  t)  of  the  relief  function,  de¬ 
fined  in  the  previous  section.  Its  Fourier  transform  obeys  the 
differential  equation 

=  -K^D  Fg(K,  t)  +  0(K.  t)  (24) 

where  12(K,  t)  is  the  spatial  Fourier  transform  of  5(R,  t).  The 
solution  is 

Fy(K,  0=  /  (25) 

Jo 

where  the  lower  integration  limit  has  been  chosen  so  that 
F^(K,  0)  =  0.  The  spectrum  for  the  growing  roughness 

follows  as 


IV,(K,  t) 

=  /  /  ti,  12)  cxp[-A"^D  {2t  -  -  t2)]dtidt2^ 

Jo  Jo 

(26) 


Here,  ^1,  ^2)  is  the  two-time  cross  spectrum  for  the 

forcing  term  5(R,  /.),  obtained  in  terms  of  the  second  moment 
of  t)  using  an  expression  analogous  to  (6).  Because 

the  forcing  term  is  stationary  with  respect  to  time,  the  cross 
spectrum  depends  only  on  r  =  ^2  —  ^1-  Changing  integration 
variables  in  (26)  to  r  and  t'  =  (ii  -f-  <2)/2,  the  integral  over  t' 
can  be  performed,  resulting  in 


W,{K,  t)  = 


K^D 


■  [  W,{K,T){cxp{-K^DT)-cxp[K‘^D{T-2f)]}dT. 

Jo 

(27) 


Dividing  the  spectrum  into  decaying  and  growing  parts  repre¬ 
senting  /a(R,  t)  and  /y(R,  t) 


ir(K.  0.  t)  =  R;.(K)  e.xp{-K^Dt.)  (22) 


U;.(K)  =  R'KK,  ()4-VK„(K,  t)  (28) 
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then 


0.02  day.  h  =  0.04  cm 


t)  =  \V„{K)exv{-2DK'^t)  (29) 

and 

W^{K.  t)  =  ^n.(K)[l  -  ^xv{-2DK^t)].  (30) 

It  can  be  seen  that  the  following  expression  for  H4(K,r)  in¬ 
serted  in  (27)  yields  (30): 


0  05  I 


y-coordinatc  (m)  x-coordinate  (m) 

0.5  day.  h  =  0  1 1  cm 


IK.(K,  r)  =  2K^D  )y„(K)^(r).  (31) 

In  this  substitution,  note  that  the  only  contribution  to  the  integral 
(27)  comes  at  r  =  0.  To  make  sense  of  the  integral  of  the  delta 
function,  it  should  be  regarded  as  having  finite  width  and  being 
an  even  function  of  r  so  that  its  integral  from  0  to  any  finite 
positive  T  yields  1/2. 

Expression  (31)  is  the  primary  result  of  this  section.  While 
equivalent  to  [5,  eq.  (5.91)],  it  has  been  derived  from  assump¬ 
tions  of  independence  and  stationarity  rather  than  from  a  spe¬ 
cific  model  for  the  forcing  term.  This  expression  shows  that  the 
two-time  cross  spectrum  of  the  forcing  term  is  completely  de¬ 
termined  by  the  naturally  occurring  roughness  spectrum  and  the 
diffusivity.  The  factor  in  (3 1 )  enhances  the  short-wavelength 
parts  of  the  forcing  spectrum  compared  to  the  natural  spectrum. 
The  delta  function  in  the  time  lag  variable  shows  that  the  forcing 
term  has  a  white  temporal  spectrum;  that  is,  successive  values 
of  the  random  forcing  time  series  are  uncorrelated.  This  is  con¬ 
sistent  with  the  statement  made  in  the  Introduction  that  the  time 
scales  associated  with  forcing  are  very  short  so  that  longer  term 
nonstationarity  could,  in  principle,  be  incorporated  by  a  rather 
straightforward  extension  of  the  present  model. 

The  forcing  spectrum  may  provide  insight  into  details  of 
the  processes  that  create  roughness.  The  factor  2K^D  is 
the  average  spectrum  of  the  features  created  per  unit  time.  In 
other  words,  if  one  considers  the  squared  magnitude  of  the  2-D 
Fourier  transform  of  the  relief  function  for  each  biologically 
created  pit,  mound,  etc.,  and  averages  these  with  weighting 
corresponding  to  rate  of  occurrence,  the  result  is  proportional 
to 

To  illustrate  the  action  of  the  forcing  term,  simulations  will  be 
used  in  which  a  random  Gaussian  forcing  time  and  space  series 
is  generated  so  as  to  obey  (31).  Equation  (24)  is  implemented 
numerically,  and  an  inverse  Fourier  transform  is  used  to  obtain 
realizations  of  the  growing  relief  function  /^(R,  t).  This  pro¬ 
vides  a  simulated  view  of  the  roughness  that  would  grow  nat¬ 
urally  if  the  seafloor  were  perfectly  flat  at  /  =  0.  The  natural 
spectrum  is  assumed  to  have  the  “von  Karman”  form 


y-coordinatc  (m)  x-coordinate  (m) 

25  days,  h  =  0.26  cm 


0.05 


Fig.  2.  Simulated  growing  naiural  relief  generated  using  the  parameters  D  = 
10“*^  m^s“\  72  =  3.5,  W2  =  0.001  m'’-’',  and  A'n  =  2  tt/O.I  m“^  The 
interface  is  perfectly  flai  at  time  zero,  and  ihe  growing  RMS  relief  is  denoted  h. 

Fig.  2  shows  three  selections  from  a  simulated  time  series  of 
roughness  realizations  generated  using  the  method  discussed 
in  the  Appendix.  The  last  of  the  three  realizations  corresponds 
to  a  time  of  25  days  after  the  initial  flat  condition.  After  this 
relatively  long  time,  roughness  is  near  equilibrium  with  RMS 
relief  0.26  cm,  to  be  compared  with  the  equilibrium  value 
0.29  cm  given  by  (33).  Close  inspection  of  Fig.  2  shows  that 
the  later  realizations  are  richer  in  longer  wavelength  roughness 
than  the  earliest  realization.  This  is  seen  more  clearly  in  the 
corresponding  spectra  presented  in  Fig.  3,  which  shows  that 
short-wavelength  parts  of  the  spectrum  grow  toward  equilib¬ 
rium  rapidly  while  the  longer  wavelength  parts  grow  more 
slowly.  In  fact,  the  factor  in  (31)  suppresses  the  growth  of 
large  features  to  such  an  extent  that  the  spectra  exhibit  a  “hole’' 
centered  at  the  origin.  This  hole  is  evident  in  the  theoretical 
spectra  of  the  left-hand  side  panels  but  is  not  obvious  in  the 
spectra  obtained  from  the  roughness  realizations,  as  the  statis¬ 
tical  uncertainty  is  rather  large. 


IV'niK)  = 


_ W2 _ 

{K^  +  ' 


(32) 


As  with  all  spectra  employed  in  this  paper,  the  integral  over  all 
K  gives  the  variance  of  the  random  variable  represented  by  the 
spectrum.  In  this  case,  the  integral  gives  the  mean  square  relief 


27r  VJ2 

(72  -  ■ 


(33) 


III.  ACOUSTIC  Measurement  Issues 

Observation  of  the  time  dependence  of  backscattering 
strength  and  ping-to-ping  correlation  involves  issues  not 
encountered  in  typical  backscatter  measurements.  Any  mea¬ 
surement  of  time  dependence  must  employ  sufficient  averaging 
to  reduce  the  effect  of  statistical  fluctuations  to  a  level  that  will 
not  mask  the  sought-after  change.  In  some  respects,  however, 
acoustic  measurements  of  roughness  evolution  are  easier  than 
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Fig.  3.  The  spectra  of  the  simulated  surfaces  of  Fig.  2  (right-hand  side  panels) 
compared  with  the  theoretical  spectra  given  by  (30)  (left-hand  side  panels). 


conventional  measurements  of  backscattering,  which  demand 
accurate  knowledge  of  experimental  geometry  as  well  as  trans¬ 
ducer  properties  including  source  level,  receiver  sensitivity, 
and  directivities.  In  contrast,  the  measurements  to  be  discussed 
in  the  following  sections  involve  comparisons  of  acoustic  data 
taken  at  different  times  and  do  not  require  complete  characteri¬ 
zation  of  the  measurement  apparatus  and  geometry. 

A  measurement  of  backscattering  strength  is  essentially  a 
measurement  of  the  mean  square  pressure  envelope.  The  re¬ 
quired  averaging  is  usually  obtained  by  means  of  successive 
pings  taken  from  different  patches  of  the  seafloor.  This  is  ac¬ 
complished  either  by  motion  of  the  sonar  over  the  seafloor  or 
by  pointing  the  sonar  successively  toward  different  portions  of 
the  seafloor.  These  methods  are  not  acceptable  in  the  present 
case  because  the  changes  in  scattering  patch  location  will  mask 
changes  due  to  biological  activity.  Instead,  the  mean  square  en¬ 
velope  can  be  estimated  from  single  echoes  by  means  of  an  av¬ 
erage  over  a  portion  of  the  echo  time  series 

I  P{t)r*{t)dt.  (34) 

Here,  P{t)  is  the  complex  (baseband)  pressure  signal,  the  as¬ 
terisk  denotes  complex  conjugation,  and  the  integral  is  over  a 
time  window  of  length  T.  This  window  is  the  “range  gate,’’  cen¬ 
tered  at  a  time  corresponding  to  the  range  of  interest.  A  nota- 
tional  difficulty  arises  at  this  point  because  two  different  time 


scales  are  required  in  this  paper.  The  time  scale  used  in  (34) 
can  be  termed  “fast  time”  as  it  is  measured  in  milliseconds  and 
is  the  round-trip  acoustic  time-of-flight  from  the  sonar  to  the 
seafloor  and  back.  This  is  distinct  from  “slow  time,”  which  is 
the  time  between  any  two  pings  that  are  to  be  compared.  This 
slow  time  was  used  in  the  model  development  earlier  and  is  of 
the  order  of  the  time  required  for  significant  changes  to  occur  in 
the  seafloor,  hours  to  days.  In  practice,  the  integral  in  (34)  will 
be  approximated  by  a  sum  over  samples  taken  at  discrete  times. 
If  the  pressure  is  a  Gaussian  random  process,  it  can  be  shown 
that  the  normalized  standard  deviation  for  this  estimate  (stan¬ 
dard  deviation  divided  by  the  mean)  can  be  approximated  by 


(4)  “  Vtb 

where  B  is  an  effective  bandwidth,  defined  as  [26] 


(35) 


immi] 


In  this  expression,  4>(/)  is  the  power  spectrum  of  the  trans¬ 
mitted  waveform.  Remembering  that  all  signals  are  represented 
in  baseband,  the  spectrum  is  centered  at  /  =  0,  and  it  can  be 
seen  that  the  effective  bandwidth  coincides  with  the  width  of  a 
rectangular  spectrum,  as  it  should.  The  approximation  leading 
to  (35)  assumes  that  the  time-bandwidth  product  TB  \s  much 
larger  than  unity.  The  symbols  (}  in  (35)  have  the  same  meaning 
as  earlier:  they  denote  a  formal  ensemble  average.  Thus,  two 
types  of  averages  are  of  interest.  Sample  averages,  as  in  (34),  are 
applied  to  data  and  result  in  estimates  of  such  quantities  as  mean 
square  pressure  and  correlation.  Formal  averages  give  the  theo¬ 
retical  predictions  of  the  model  and  are  also  used  in  estimating 
the  error  in  sample  averages,  as  in  (35).  The  time-bandwidth 
product  can  be  interpreted  as  the  effective  number  of  indepen¬ 
dent  samples  in  the  windowed  time  series.  Thus,  the  uncertainty 
in  estimates  such  as  (34)  is  not  determined  by  the  actual  number 
of  samples,  as  an  increase  in  sampling  rate  results  in  higher  cor¬ 
relation  between  adjacent  samples  and  does  not  improve  the  es¬ 
timate.  In  consequence,  the  sampling  rate  is  not  of  particular 
importance  as  long  as  the  Nyquist  criterion  is  satisfied. 

If  one  hopes  to  observe  changes  in  scattering  strength  on  the 
order  of  3  dB,  the  normalized  standard  deviation  must  be  smaller 
than  unity.  Thus,  it  is  desirable  to  have  large  bandwidth,  though 
not  so  large  that  the  single-frequency  assumption  inherent  in  the 
above  perturbation-theory  expressions  is  not  satisfied.  The  aver¬ 
aging  window  width,  T  sets  the  range  resolution,  and  it  follows 
that  one  is  better  off  not  seeking  high  resolution  by  reduction  in 
the  averaging  window  width.  The  value  of  1/(TZ?)*/^  for  the 
acoustic  results  to  be  reported  here  is  about  0.27  for  processing 
of  the  40-kHz  data  and  0.39  for  processing  of  the  300-kHz  data. 

Even  for  relatively  “hard”  sandy  seafloors,  backscattering  is 
not  wholly  due  to  interface  roughness;  there  will  be  a  volume 
scattering  component  resulting  from  heterogeneity  of  the  sed¬ 
iment.  This  situation  can  be  accommodated  in  (19)  by  identi¬ 
fying  (Tn{0)  as  the  natural  scattering  cross  section  resulting  from 
both  roughness  and  sediment  heterogeneity.  In  short,  (7„(^)  is 
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the  asymptotic  value  of  the  cross  section  that  is  observed  be¬ 
fore  the  artificial  roughness  is  created  and  after  it  has  decayed 
to  zero. 

Correlation  estimates  can  be  made  using  the  following  com¬ 
plex  correlation  coefficient  [12],  [14],  [15],  [27],  [28]: 


C(r)  = 


f  P2(t.)P{(t)dt. 


lfP,(f.)P,^(t)dtfP,(t)Pi(t)dt] 


1/2 


(37) 


where  Pi(t)  and  P2(t)  are  the  complex  baseband  echo  signals 
received  from  two  different  pings  transmitted  with  time  differ¬ 
ence  T  =  t'z  —  ti.  where  these  are  “slow  times.”  The  integral 
is  over  a  (fast  time)  window  of  width  T.  The  correlation  esti¬ 
mate  provided  by  (37)  is  an  irrational  function  of  the  signals 
and  rather  difficult  to  analyze  statistically.  In  order  to  approxi¬ 
mate  the  statistical  “floor”  of  the  correlation  estimate,  one  can 
consider  the  case  of  two  echo  signals  P\{t)  and  P2(t)  that  are 
Gaussian  and  uncorrelated,  with  TB  ^  1.  In  this  case,  the  ap¬ 
proximate  formal  expression 


Kf  ~tb 


(38) 


can  be  obtained,  where 

Ki2  =  ^  I  P2{t)Pr{t)dt.  (39) 

Expression  (38)  is  not  the  standard  deviation  of  the  correlation 
estimate  (for  two  uncorrelated  signals),  because  the  required 
variance  would  have  to  deal  with  moments  of  the  irrational  ex¬ 
pression  (37).  However,  if  TB  ^  1,  the  relative  fluctuations 
of  the  normalizing  denominator  of  (37)  will  be  small.  If  these 
fluctuations  are  neglected,  (38)  can  be  considered  an  approxi¬ 
mation  to  the  variance  of  the  correlation  estimate  for  the  case 
of  two  uncorrelated  echoes.  If  one  hopes  to  estimate  nonzero 
correlations  that  are  substantially  less  than  unity,  the  standard 
deviation  of  the  estimate  must  be  small  compared  to  unity,  that 
is,  l/{TBy^‘^  <C  1.  In  the  work  to  be  reported,  correlation  esti¬ 
mates  are  made  for  each  pixel  in  a  sonar  image  (with  range-gate 
width  T)  using  the  absolute  value  of  (37).  These  absolute  values 
are  then  averaged  over  many  pixels  to  provide  the  final  correla¬ 
tion  estimate.  Even  if  the  echoes  being  compared  are  completely 
uncorrelated,  the  estimate  will  be  nonzero,  and  this  “floor”  value 
can  be  approximated  as  follows.  First,  as  the  estimate  (37)  in¬ 
volves  a  sum  with  several  terms,  it  will  be  approximated  as 
Gaussian,  with  absolute  value  that  is  Rayleigh-distributed.  Ac¬ 
cording  to  (38),  the  mean  square  value  of  this  Rayleigh  random 
variable  is  approximately  l/(TB).  It  follows  that  the  asymp¬ 
totic  value  of  the  mean  correlation  estimate  is  approximately 

(|C(rxD)|>  =  i\/^.  (40) 


Just  as  in  the  case  of  scattering  strength  estimates,  it  is  desirable 
to  have  large  bandwidth  and  as  wide  an  averaging  time  window 
as  permitted  by  range  resolution  requirements. 

As  with  scattering  strength  estimates,  correlation  estimates 
are  complicated  by  the  fact  that  both  roughness  and  volume 
scattering  contribute  to  the  echo  signal.  Because  volume  echoes 


tend  to  decorrelate  on  much  longer  time  scales  than  those  due  to 
roughness  scattering  [18],  a  simple  expression  that  can  accom¬ 
modate  volume  scattering  can  be  obtained  by  a  slight  general¬ 
ization  of  (23)  and  (37)  leading  to 

C(t)  =  (l-7)cxi)  +7  (41) 

where  0  <  7  <  1  is  the  asymptotic  value  of  the  correlation  es¬ 
timate  for  large  lag  times,  and  where  Ta  =  2T£),  with  Tp  given 
by  (16).  The  nonzero  value  of  this  asymptote  is  due  to  both  the 
presence  of  unchanging  scattering  (e.g.,  due  to  volume  hetero¬ 
geneity  not  subject  to  rapid  change)  and  also  due  to  the  unavoid¬ 
able  variance  of  the  correlation  estimate  as  discussed  above. 
Thus,  the  parameter  7  must  be  larger  than  (1(7(00)1)  as  given 
by  (40).  Expression  (41)  is  presented  as  a  hypothesis  without 
proof,  and  the  parameter  7  has  the  burden  of  incorporating  both 
unknown  scattering  mechanisms  and  the  statistical  “floor”  of 
correlation  estimates  (40).  This  parameter  will  be  determined 
by  fits  to  ping-to-ping  correlation  data. 

IV.  ACOUSTIC  Measurement  Systems 

Two  bottom-mounted,  rotating  sonar  systems  were  used, 
the  BAMS,  operating  at  40  kHz  and  the  XBAMS,  operating 
at  300  kHz.  These  are  autonomous  systems  that  execute  cir¬ 
cular  scans  according  to  a  preset  schedule.  The  azimuthal, 
3-dB,  full  beamwidths  were  5°  for  the  40-kHz  system  and 
1°  for  the  300-kHz  system.  These  beamwidths  are  equal  to 
the  azimuthal  step  made  between  successive  pings.  As  the 
same  transducer  was  used  for  transmission  and  reception,  the 
round-trip  beamwidths  are  smaller  by  a  factor  of  about  0.71, 
and  the  sidelobe  levels  are  quite  low.  Both  systems  transmitted 
unshaded  linear  frequency  modulated  up-sweeps  having  source 
levels  of  approximately  217  dB  re  1  //Pa  and  pulse  lengths  of 
2  ms.  The  sweep  widths  were  2  and  10  kHz  centered  at  40  and 
300  kHz,  respectively.  The  echo  data  analyzed  in  this  paper 
were  in  complex  baseband  form  with  sampling  frequencies 
of  2  and  20  kHz.  These  data  were  matched  filtered  before 
final  processing  with  the  exception  of  the  earliest  experiment 
[Coastal  Benthic  Boundary  Layer  (CBBL),  reported  below]. 
Matched  filtering  has  little  effect  on  the  40-kHz  cross-section 
and  correlation  estimates,  as  the  window  width  (T  =  6.7  ms)  is 
large  compared  to  the  pulse  length  (2  ms). 

Both  BAMS  and  XBAMS  employed  tripod  platforms.  For 
the  earlier  BAMS  measurements  reported  below  (CBBL),  the 
height  of  the  sonar  above  the  seafloor  was  5.2  m.  For  the  later 
SAX99  and  SAX04  measurements,  the  transducer  height  for 
both  BAMS  and  XBAMS  was  3  m. 

V.  Measurement  Sites 

The  data  to  be  compared  with  the  diffusion  model  come  from 
the  CBBL  Panama  City  experiment  as  well  as  from  the  SAX99 
and  the  SAX04.  These  data  represent  a  small  part  of  the  total 
acoustic  and  photographic  data  from  these  experiments.  All 
three  experimental  sites  were  located  on  inner  shelf  sands  in 
the  northeastern  Gulf  of  Mexico  in  water  depths  of  15-29  m. 
The  primary  differences  in  the  sediments  at  these  three  sites 
were  in  the  percentage  of  shell  fragments  (percent  carbonate). 
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presence  or  absence  of  mud  lenses,  and  oceanographic  con¬ 
ditions  that  can  alter  seafloor  roughness.  Sediment  physical 
properties  (mean  grain  size,  porosity,  and  bulk  density)  were 
determined  from  sediment  samples  collected  using  diver  cores. 
Sediment  geoacoustic  properties  (sound  speed  and  attenuation) 
were  either  measured  in  situ  at  38,  40,  or  58  kHz,  near  the 
40-kHz  operational  frequency  of  BAMS,  or  at  400  kHz  on 
diver-collected  core  samples,  near  the  300-kHz  operational 
frequency  of  XBAMS. 

During  the  1993  CBBL  Panama  City  experiment,  the  BAMS 
sonar  platform  was  deployed  at  29°  41.40'  N,  85°  40.71'  W 
in  water  depth  of  29  m  [18].  Acoustic  scans  were  acquired  at 
20-min  intervals  over  an  18.3-day  period  (August  13-31,  1993). 
Based  on  sidescan  sonar  imagery  and  surface  grab  samples,  sed¬ 
iments  at  this  site  varied  from  highly  reflective  coarse  sand  to 
less  reflective  fine  sand  [29],  [30].  BAMS  was  located  in  an  area 
of  high  relative  acoustic  reflectivity,  with  the  sediment  com¬ 
posed  of  a  mixture  of  medium  to  coarse  sand  and  carbonate 
shell  hash  (carbonate  percentage  21%  ±12.5%).  Analysis  of 
diver-core  samples  collected  within  the  acoustic  field  of  view 
gave  a  mean  sediment  grain  diameter  of  0.48  mm  ( 1 .05  0),  mean 
porosity  41.2%,  and  mean  bulk  density  of  2000  kg  m“^.  These 
coarse-grain  sediments  had  a  mean  sound-speed  ratio  of  1.132 
and  a  mean  attenuation  of  33.7  dB  m~^,  measured  in  situ  at 
58  kHz  [30].  No  strong  bottom  currents  or  storms  with  high 
significant  wave  heights  were  observed  during  the  August  1993 
deployment  of  BAMS,  suggesting  any  changes  in  bottom  rough¬ 
ness  during  the  acoustic  experiment  were  a  result  of  biological 
activity. 

During  SAX99,  BAMS  was  deployed  approximately  150  m 
west  of  the  main  tending  vessel  at  30°  22.66'  N,  86°  38.79'  W 
in  water  depth  of  18.5  m  [23].  Acoustic  scans  were  acquired 
at  90-min  intervals  from  October  6  to  22,  1999  and  at  30-min 
intervals  over  a  15-day  period  beginning  October  22,  1999. 
Only  data  from  the  second  period  are  considered  here.  During 
the  second  period,  a  variety  of  manipulative  experiments 
were  conducted  within  the  acoustic  field  of  view  of  BAMS, 
including  placing  discrete  scatterers  on  the  seafloor  and  cre¬ 
ating  artificial  roughness  by  raking  the  seafloor  [6],  [22], 
[23],  [31].  Stereophotographs  which  were  used  to  measure 
temporal  changes  in  roughness  after  raking  were  made  near  the 
tending  vessel  November  4-6,  1999.  Based  on  backscattering 
from  3()0-kHz  multibeani  sonar  images,  the  site  around  the 
BAMS  tower  was  much  more  uniform  than  during  the  1993 
experiments.  From  analysis  of  diver-collected  core  samples, 
the  sediments  were  medium  sands  with  a  mean  diameter 
of  0.42  mm  (1.27  0),  mean  porosity  of  36.6%,  and  mean 
bulk  density  of  2074  kg  m“^  [22].  These  sediments  had  a 
mean  sound-speed  ratio  of  1.16  [31]  and  mean  attenuation  of 
12.7  dB  m“^,  measured  in  situ  at  38  kHz  [22].  The  acoustic 
measurements  reported  in  this  paper  for  SAX99  were  collected 
during  periods  of  relative  calm  (October  23-24  and  October 
30-31)  where  bottom  currents  and  significant  wave  heights  of 
surface  gravity  waves  were  well  below  the  threshold  of  motion 
for  medium-sand-sized  particles  found  on  the  seafloor.  The 
bottom  stereophotographs  reported  in  Section  VI  were  also 
collected  during  periods  of  relative  calm  with  bottom  currents 
less  than  10  cm  s“^  and  significant  wave  heights  less  than  1  m. 


These  data  suggest  that  any  changes  in  bottom  roughness  during 
the  acoustic  experiments  or  bottom  roughness  characterization 
were  the  result  of  biological  modification  of  the  seafloor. 

Experiments  were  conducted  during  SAX04  at  a  site  approx¬ 
imately  2  km  north  of  the  SAX99  experimental  site  in  water 
depth  of  16.5  m.  BAMS  was  deployed  on  September  19,  2004 
at  30°  23.284'  N,  86°  38.496'  W  and  continued  to  collect  data 
at  1-h  intervals  until  October  6,  2004.  The  deployment  was 
about  ten  days  after  landfall  of  Hurricane  Ivan,  and  the  seafloor 
was  still  covered  with  mud,  presumably  derived  from  the  back¬ 
wash  of  lagoonal  sediments.  Resuspension  of  this  mud  during 
the  11 -day  acoustic  deployment  resulted  in  poor  visibility, 
restricting  diver-manipulative  experiments  and  photography; 
however,  the  mud  probably  insulated  the  rippled  sand  seafloor 
from  the  effects  of  fish  feeding.  After  an  unsuccessful  first 
deployment,  XBAMS  was  redeployed  on  October  22,  2004  at 
30°  23.213'  N,  86°  38.782'  W  and  collected  data  until  October 
25,  2004  at  1-h  intervals.  Based  on  sidescan  imagery,  the  sites 
around  both  BAMS  and  XBAMS  towers  were  much  less  uni¬ 
form  than  during  the  1993  and  1999  experiments.  Areas  of  low 
backscatter  (presumably  mud)  were  interspersed  with  strong 
reflections  from  rippled  seafloor.  Sediments  contained  layers  of 
sand  and  layers  of  mud  derived  from  lagoonal  sediments  [33]. 
Near  the  BAMS  tower,  the  lagoonal  mud  was  draped  over  a 
rippled  sand  seafloor.  The  XBAMS  data  should  not  be  greatly 
affected  by  mud  as  they  were  acquired  after  the  mud  layer  had 
dissipated.  The  medium-sand  sediments  had  a  mean  diameter 
of  0.36  mm  (1.46  0),  a  mean  porosity  of  36.6%,  and  a  mean 
bulk  density  of  2064  kg  m”"^  [24],  [34].  These  sediments  had 
a  mean  sound-speed  ratio  of  1.161  and  a  mean  attenuation  of 
10  dB  m“^,  both  measured  in  situ  at  40  kHz  [35].  Sound-speed 
ratio  and  attenuation  measured  at  400  kHz  from  diver-core 
samples  were  1.162  and  92  dB  m“\  respectively.  The  mud 
sediments  that  covered  the  rippled  seafloor  had  a  mean  grain 
size  of  10  0,  porosity  of  80%-90%  and  a  mean  bulk  density  of 
1 150-1250  kg  m“^.  Sound  speeds  measured  at  400  kHz  were 
slightly  less  than  the  overlying  water. 

VI.  Artificial  Roughness  Measurements 

Backscattering  measurements  on  an  artificially  roughened 
seafloor  were  carried  out  as  part  of  the  SAX99  experiment  [6]. 
In  one  set  of  measurements  the  sediment  was  raked  to  produce 
regular  ripples  of  approximate  wavelength  1.95  cm,  a  value 
chosen  to  correspond  to  the  Bragg  wavelength  at  40  kHz.  As 
will  be  seen,  these  ripples  caused  elevated  levels  of  backscat¬ 
tering,  and  it  is  the  time  decay  of  this  elevated  backscattering 
that  is  of  interest.  Fig.  4  shows  images  taken  from  stereopho¬ 
tographs  of  the  seafloor  before  and  after  raking.  A  steady 
decay  of  the  artificial  roughness  is  evident,  but,  contrary  to 
the  picture  presented  in  the  Introduction,  some  of  this  decay 
resulted  from  the  activity  of  larger  bottom-dwelling  fauna  such 
as  starfish,  crabs,  sand  dollars,  and  fish,  of  which  the  flounder 
seen  in  the  last  picture  in  the  series  is  representative.  Roughness 
degradation  of  this  sort  may  not  be  diffusive  on  the  small  scales 
assumed  in  the  present  model.  The  roughness  spectra  obtained 
from  these  stereophotographs  are  shown  in  Fig.  5.  Note  that 
spatial  frequency  is  used  as  the  spectral  variable.  It  is  equal  to 
the  wave  number  used  in  other  parts  of  this  paper  divided  by 
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(C)  (d) 


Fig.  4.  Photographs  taken  91  cm  from  the  seafloor  during  SAX99  using  a  Kodak  DC  120  digital  camera  with  a  resolution  of  1280  x  960  pixels.  The  photographs 
show  the  seafloor  under  the  camera  frame  (a)  before  manipulation,  (b)  6  min  after  raking, l.95-cm  ripple  wavelength,  (c)  3  h  after  raking,  and  (d)  7  h  after  raking. 
A  flounder,  partially  responsible  for  the  decay  of  the  artificial  changes  in  the  seafloor,  can  be  seen  in  the  lower  right-hand  side  corner  of  the  7-h  photograph.  These 
photographs  were  taken  during  a  period  of  calm  weather.  (Reproduced  from  [6].) 


2tt.  Spectral  peaks  due  to  the  artificial  ripple  are  clearly  visible 
in  the  spectrum  immediately  after  raking,  but  these  are  less 
evident  after  3  h  and  not  seen  after  7  h. 

As  an  application  of  the  artificial  roughness  model  developed 
earlier,  a  value  for  diffusivity  will  be  obtained  from  the  series  of 
spectra  shown  in  Fig.  5.  The  spectrum  immediately  after  raking 
will  be  modeled  in  a  generalization  of  the  form  used  in  [6].  To  in¬ 
sure  that  the  spectrum  has  the  proper  symmetry,  it  will  be  broken 
into  parts  with  peaks  at  positive  and  negative  Ky 

\V{K)  =  W^{K)  +  W_{K) .  (42) 


These  two  parts  will  each  be  divided  into  three  subcomponents 


U/±(K)  =  ^ 

71-1 


2 

alKl  +  (Ky  ±  K,,f  + 


(43) 

The  parameters  Kn  allow  introduction  of  spectral  peaks,  and 
the  parameters  allow  control  of  the  anisotropy  between  the  x- 
and  ?y-directions  in  wave  vector  space.  The  first  term  in  (43)  will 
be  taken  to  represent  the  natural  background  isotropic  rough¬ 
ness,  with  the  parameters  ATi,  ai,  /i,  W2iy  and  721  set  to  0,  1, 
3  cm,  0.05  cm”\  and  5,  to  provide  a  rough  match  to  the  spec¬ 
trum  before  raking.  The  second  term  in  (43)  represents  the  arti¬ 
ficial  roughness  at  low  spatial  frequencies  with  iir2>  a2,  /2,  ^22^ 
and  722  set  to,  0,  6,  3  cm,  0.25  cm“^,  and  5.  The  third  term 


has  A'3,  03,  /3,  W2'3,  and  723  set  to  27r/1.95  cm“^,  0.7,  7  cm, 
0.05  cm^*^,  and  2.5  to  match  the  spectral  peak  due  to  raking.  To 
employ  the  model  developed  for  decay  of  artificial  roughness, 
it  is  only  necessary  to  append  the  factor  exp{~2K^ Dr)  to  the 
last  two  terms  in  (43),  as  these  are  the  nonequilibrium  terms 
subject  to  decay.  The  resulting  time-dependent  model  spectrum 
is  shown  in  Fig.  6  for  the  times  given  in  Fig.  4.  The  diffusivity 
parameter  was  taken  to  be  D  —  3  x  10“^^  m^  s“^  in  order  to 
obtain  an  approximate  match  with  Fig.  5.  A  change  in  the  as¬ 
sumed  diffusivity  by  25%  produced  a  clear  mismatch.  Fig.  6 
illustrates  an  important  feature  noted  earlier:  high-frequency 
(short-wavelength)  features  decay  more  rapidly  than  low-fre¬ 
quency  features.  Thus,  the  two  spectral  peaks  near  50  cycles/m 
due  to  raking  disappear  before  the  low-frequency  lobe  has  re¬ 
turned  to  its  natural,  isotropic  form. 

During  SAX99,  backscattering  measurements  [6]  were  made 
using  BAMS  on  artificial  ripple  fields  similar  to  that  charac¬ 
terized  by  the  spectra  of  Fig.  5.  These  ripples  were  oriented 
with  crest  lines  (equivalently,  “strike”)  normal  to  the  direction 
of  the  incident  sound.  The  ripple  wavelength  of  1.95  cm  was 
very  close  to  the  wavelength  for  Bragg  scattering,  1.99  cm, 
which  results  for  the  frequency  (40  kHz),  grazing  angle  (15°), 
and  water  sound  speed  (1538  m  s“^)  of  this  measurement. 
Scattering  strength  was  obtained  by  averaging  the  scattering 
cross  section  over  all  pixels  within  the  2-m  x  2-m  raked  area, 
with  center  at  a  horizontal  range  of  1 1  m  from  the  BAMS  rota¬ 
tion  axis.  This  spatial  average  was  then  converted  to  scattering 
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Fig.  5.  Two-dimensional  height  speciral  density  levels  (dB  re  m'*)  calculated  from  ihe  siereo-pair  phoiographs  corresponding  to  the  images  depicted  in  Fig.  4. 
The  peaks  near  50  cycles/m  are  due  lo  raking  of  ihe  seafloor.  (Reproduced  from  [6].) 


Strength,  which  is  equal  to  101ogiQfr(6/)  where  a{0)  is  the 
scattering  cross  section  defined  in  (18).  Fig.  7  shows  a  model 
fit  to  the  results  of  two  separate  rakings,  with  abrupt  increase  in 
scattering  strength  and  subsequent  decay.  The  two  rakings  have 
differing  scattering  strengths,  but  the  decay  rate  of  both  is  sim¬ 
ilar.  The  model  fit  results  in  £)  =  1.8  x  m^  s“^  Fig.  7 

shows  the  result  of  25%  changes  in  diffusivity,  cited  earlier  as 
the  uncertainty  in  fitting  the  spectra  of  Fig.  5.  The  diffusivity 
obtained  by  fitting  backscattering  data  for  SAX99  is  a  factor  of 
1 7  smaller  than  the  value  obtained  by  fitting  the  decaying  spec¬ 
trum.  The  difference  in  diffusivity  values  obtained  from  decay 
of  the  relief  spectrum  versus  decay  of  backscattering  accords 
with  qualitative  observations.  The  relief  spectrum  was  obtained 
from  a  site  near  the  SAX99  tending  vessel,  which  acted  as  an 
artificial  reef  and  attracted  numerous  fish  (Fig.  1).  The  BAMS 
acoustic  measurements  were  conducted  a  few  hundred  meters 
from  the  vessel  where  feeding  fish  were  not  so  numerous.  If 
the  feeding  activity  accounts  for  roughness  decay  as  well  as 
roughness  creation,  this  could  explain  the  large  difference  in 
diffusivity  seen  at  the  two  sites.  This  comparison  makes  it  clear 


that  future  experiments  should  employ  roughness  and  acoustic 
measurements  that  are  coincident  in  both  position  and  time.  If 
the  resulting  diffusivities  agree  over  a  wide  range  of  forcing 
conditions,  this  would  be  a  partial  validation  of  the  model. 

VII.  Natural  Roughness  Measurements 

Ping-to-ping  correlation  measurements  on  natural  sandy 
seafloors  were  carried  out  during  three  experiments  in  the 
northeastern  Gulf  of  Mexico,  the  CBBL  experiment,  SAX99, 
and  SAX04.  During  all  three  experiments,  measurements  were 
made  at  40  kHz  using  BAMS.  During  SAX04,  measurements 
were  also  made  at  300  kHz  using  XBAMS.  Correlation  es¬ 
timates  were  made  for  adjacent,  nonoverlapping  averaging 
windows  and  for  each  azimuthal  position  in  the  circular  scans 
provided  by  BAMS  and  XBAMS.  Each  estimate  can  be  thought 
of  as  providing  a  value  for  a  pixel  in  a  correlation  image  of  the 
seafloor  near  the  sonar.  Finally,  the  correlation  (absolute)  values 
were  averaged  over  a  large  number  of  adjacent  pixels,  taking 
care  to  avoid  inclusion  of  regions  of  obviously  different  scat¬ 
tering  strength  or  correlation  decay.  The  results  reported  here 
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X-frequency  (m"  )  X-frequency  (m’  ^ ) 


Fig.  6.  Model  spectra  for  comparison  with  the  measured  spectra  of  Fig.  5.  The  assumed  diffusivity  is  O  =  3  x  10“^  m^ 


Fig.  7.  Model  fit  to  decay  of  scattering  strength  alter  two  different  raking  treat¬ 
ments  during  SAX99.  The  sonar  frequency  is  40  kHz  and  the  nominal  grazing 
angle  is  15°.  The  solid  curve  has  D  =  1.8  x  10“'*^  s“ ‘ .  and  the  dashed 
curves  correspond  to  diffusivities  larger  by  a  factor  of  1 .25  and  smaller  by  a 
factor  of  0.75.  The  initial  scattering  strength  is  assumed  to  be  —24  dB  and  the 
asymptotic  scattering  strength  is  assumed  to  be  —  36  dB. 


are  merely  a  sampling  from  the  available  circular  correlation 
images.  For  the  40-kHz  measurements,  the  bandwidth  D  was 
2  kHz  and  the  averaging  window  width  T  was  chosen  such  that 
(T /2  =5.0  m.  For  the  300-kHz  measurements,  B  =  \0  kHz, 
and  cT/2  =  0.5  m.  The  corresponding  asymptotic  correlation 
estimates  (39)  are  0.24  for  BAMS  and  0.34  for  XBAMS. 
These  figures  represent  the  statistical  “floor'  for  correlation 
estimates  for  large  lag  times.  As  noted  above,  the  actual  floor  7 
will  be  higher  if  slowly  changing  sediment  volume  scattering 
contributes  to  the  echo  signal. 


Fig.  8  shows  40-kHz  correlation  data  and  model  fits  for 
CBBL  and  SAX04.  The  CBBL  data  were  processed  using 
the  same  averaging  window  width  as  for  SAX99  and  SAX04. 
The  CBBL  correlation  estimates  were  averages  over  all  72 
azimuthal  steps  at  a  horizontal  range  of  27.5  m.  The  SAX99 
and  SAX04  correlation  estimates  were  averages  over  approx¬ 
imately  30  azimuthal  steps  centered  at  a  horizontal  range  of 
12.5  m.  The  fitted  asymptotic  correlation  (7)  was  0.4  for  CBBL 
and  0.45  for  SAX04.  These  asymptotes  are  larger  than  the 
statistical  floor  of  0.24;  consequently,  it  appears  that  roughly 
7%-20%  of  the  backscattered  intensity  must  be  due  to  sediment 
volume  scattering  that  undergoes  little  change  on  time  scales 
of  several  days.  The  SAX04  correlation  curve  is  rather  uneven 
owing  to  the  intrusion  of  fish  into  the  sonar  field  of  view.  This 
problem  was  even  more  severe  for  SAX99,  as  evident  in  Fig.  9. 
In  this  case,  the  diffusivity  given  in  Section  VI  for  scattering 
strength  decay  at  this  site  provides  a  satisfactory  fit.  Fig.  10 
compares  measured  and  modeled  ping-to-ping  correlation  at 
300  kHz.  These  correlation  estimates  were  averages  over  ten 
azimuthal  steps  and  four  horizontal  range  increments  (0.5  m 
in  size)  centered  at  a  horizontal  range  of  12  m.  The  decay 
rate  of  correlation  varied  substantially  over  the  circular  field 
of  view  of  XBAMS,  and  the  resulting  diffusivities  vary  over 
roughly  an  order  of  magnitude,  being  largest  near  the  NURC 
camera  system,  which  attracted  fish.  The  asymptotic  values  of 
the  correlation  estimates  (0.45  and  0.5)  are  somewhat  larger 
than  the  statistical  floor  (0.34),  suggesting  a  small  scattering 
contribution  that  is  relatively  constant  with  time.  There  is 
some  question  as  to  whether  roughness  or  volume  scattering 
dominates  at  300  kHz  [36].  Volume  scattering  presumably 
occurs  very  near  the  interface,  as  the  incident  grazing  angle 
is  well  below  the  critical  angle.  If  the  interface  were  perfectly 
flat,  the  evanescent  pressure  field  in  the  sediment  would  have 
an  e-folding  depth  of  only  2  mm.  Thus,  it  may  be  that  volume 
scattering  features  at  such  shallow  depths  are  subject  to  rapid 
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Fig.  8.  Comparison  of  measured  and  modeled  ping-io-ping  correlation  at 
40  kHz  for  two  sandy  sites  off  the  Rorida  Panhandle:  CBBL  and  SAX04. 
The  smooth  curves  are  model  results  for  CBBL  (D  =  1.3  x  s^V 

=  0.4.  f)  =  10.7°)  and  SAX04  {D  =  9  X  10"“  s-^  7  =  0.45, 

a  =  13.5°). 


^^0  1  2  3  4  5  6  7 


Lag  (hours) 

Fig.  10.  Comparison  of  measured  and  modeled  ping-to-ping  correlation  at 
300  kHz  for  SAX04.  The  three  smooth  curves  are  model  results  for  D  =  3.5 
X  10”^’  m^  s’’ ,  7  =  0.5  (upper  curve);D  =  6  x  10“ ”  m^  s“’,  7  =  0.5 
(middle  curve);  and  D  =  2.5  x  10“’°  m^  s“\  7  =  0.45  (lower  curve).  All 
curves  are  computed  for  a  grazing  angle  of  14°. 


Lag  (days) 


Fig.  9.  Comparison  of  measured  and  modeled  ping-to-ping  correlation  at 
40  kHz  for  the  sandy  SAX99  site.  The  smooth  curve  is  a  model  result  for  D  = 
1.8  X  10'’^  m'^  s-’,-  =  0.45,^  =  13.5°. 

change.  Even  so,  it  seems  plausible  that  this  change  could  be 
described  as  a  diffusive  process. 

The  SAX04  and  SAX99  sites  are  in  close  proximity  so  it  is 
reasonable  to  compare  the  diffusivities  obtained  in  these  two 
experiments.  For  SAX99,  the  diffusivity  obtained  from  either 
backscattering  decay  due  to  artificial  ripples  or  decorrelation 
of  natural  roughness  was  1.8  x  10~^®  m^  s“^  For  SAX04, 
diffusivities  are  somewhat  smaller,  9  x  10”^^  m^  s“^  at  40  kHz 
and  about  5  x  10“^^  m^  s”^  at  300  kHz  at  locations  where 
fish  populations  were  relatively  small.  It  is  encouraging  that  the 
SAX04  diffusivities  measured  at  40  and  300  kHz  are  of  similar 
magnitude  even  though  the  squared  frequencies  (which  appear 
in  the  inversion  equation)  are  different  by  a  factor  of  56.  There 
is  no  reason  to  expect  these  diffusivities  to  be  identical,  as  the 


measurements  were  made  at  separate  locations  and  at  different 
times.  The  diffusivities  measured  by  Hay  [2]  using  observations 
of  the  decay  of  naturally  generated  ripples  in  SAX04  are  larger, 
being  in  the  range  of  1-20  x  10”^  m^  s”^.  As  noted  by  Hay  [1  ], 
his  apparatus  attracted  numerous  fish,  resulting  in  rapid  ripple 
decay. 

VIII.  Discussion 

Comparisons  between  diffusivities  measured  at  different 
frequencies  or  by  different  methods  were  complicated  by  three 
factors.  First,  the  tending  vessel  attracted  numerous  fish,  so  that 
sites  near  the  ship  had  more  intense  biological  activity  than 
those  at  a  greater  distance.  On  a  smaller  scale,  sonar  and  camera 
tripods  also  attracted  fish,  though  the  BAMS  measurements 
were  made  at  distances  of  order  10  m  from  the  BAMS  tripod, 
where  fish  were  relatively  sparse.  Finally,  raking  the  sediment 
exposed  benthic  organisms  and  attracted  feeding  fish.  The  least 
bias  owing  to  disturbance  of  the  natural  environment  should 
result  if  correlation  measurements  are  made  at  ranges  of  tens  of 
meters  from  an  autonomous  platform,  well  removed  from  any 
tending  vessels. 

The  diffusivity  values  obtained  in  this  work  are  reasonable, 
but  further  investigation  is  needed  to  validate  or  improve  the 
model.  Future  experiments  should  include  ground-truth  mea¬ 
surement  of  diffusivity,  e.g.,  by  tracers,  photography,  or  high- 
resolution  rotary  imaging  sonar  with  short-term  deployment  to 
avoid  attraction  of  fish.  Also,  multiple  frequency  measurements 
should  be  made  at  the  same  time  to  see  if  the  estimated  diffu¬ 
sivity  is  frequency  independent  as  the  model  assumes.  It  is  pos¬ 
sible  to  exclude  larger  fauna  to  eliminate  much  of  the  forcing 
term.  If  the  forcing  term  is  eliminated,  the  diffusive  term  re¬ 
mains  and  roughness  should  decay,  at  least  at  frequencies  of  in¬ 
terest.  More  detailed  modeling  of  the  forcing  term,  as  in  [5], 
could  be  used  to  incorporate  knowledge  of  the  size  and  shape 
of  feeding  pits,  etc.  Models  for  volume  transport  such  as  those 
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of  [37]  and  [38]  could  conceivably  be  adapted  to  the  problem 
of  roughness  evolution  and  could  be  used  to  perform  numer¬ 
ical  experiments  suitable  for  testing  model  assumptions  such  as: 
1 )  the  assumed  identity  between  volume  horizontal  diffusivity 
and  roughness  diffusivity,  and  2)  independence  of  the  genera¬ 
tion  and  decay  mechanisms. 

Appendix 

The  simulation  results  displayed  in  Figs.  2  and  3  were  ob¬ 
tained  by  means  of  a  fast  Fourier  transform  (FFT)  implemen¬ 
tation  of  (24).  The  evolving  relief  function  is  represented  in 
the  spatial  frequency  domain  by  its  discrete  Fourier  transform 
F,{7n,  7i)  with  the  spatial  frequency  indices  m  and  n  ranging 
from  1  to  iV,  with  zero  frequency  corresponding  to  m  =  n  =  ly 
and  with  spatial  wave  number  sampling  interval  AK  =  27r/L, 
where  the  surface  to  be  simulated  has  width  L  in  both  x-  and 
^-dimensions.  The  time  index  is  denoted  i  with  time  sampling 
interval  The  equation  for  time  evolution  (24)  takes  on  the 
discrete  form 


Fi^i{7n,  7i)  =  Fi(m,  7i)  exp{-DKf^^^At)  +  ^^(m,  7i)At 

(Al) 

with 


7l)  =  a^nnii) 


2DI<ljVr,{I<^n.  Kyn) 


-,1/2 


(AK^At) 


.  (A2) 


The  coefficients  (/)  are  complex  Gaussian  random  variables 
with  the  conjugate  symmetry  required  to  insure  that  the  inverse 
FFT  of  Fe(m.  /i)  is  real.  To  produce  the  desired  white  temporal 
spectrum,  (i„ui('i)  and  aj,y(j)  are  independent  if  i  ^  j.  For  the 
iih  time  step,  independent,  complex  random  numbers 
are  generated,  with  the  real  and  imaginary  parts  of  being 

independent  and  normally  distributed.  The  a„in(0  formed 
from  the  as  follows: 

~  ^(A/'-m-j-2)(.V-n  +  2)(0]  *  (A3) 

The  indexing  convention  in  the  conjugated  term  assumes  period¬ 
icity  in  each  index  with  period  N.  For  example,  setting  m  =  1, 
N  -  711  2  =  N  +  I  =  1,  modulo  N.  Thus,  the  1 1  element 

of  is  equal  to  the  real  part  of  In  the  same  way, 

one  finds  that,  for  rn  =  1  and  77.  >  1,  -1-2) (.v -71+2)  = 

/q(,v-7i4-2)*  This  says  that  the  first  row  of  is  to  be  flipped 

in  the  left-right  sense,  with  the  second  column  replaced  by  the 
iVth,  the  iVth  by  the  second,  and  so  forth.  The  first  column  is 
treated  analogously,  and  the  remaining  (A”  —  l)x(7V’  —  1)  ma¬ 
trix  is  flipped  in  both  indices. 

Using  (A2)  and  (A3),  and  the  normality  and  independence  of 
the  random  numbers,  the  discrete  spatial  spectrum  of  the  forcing 
term  can  be  obtained  from 


(^2i(m.  7i)^l*{7ri\  71^)) 

=  2D{Kl,  + 


^771777/  ^71  n' 

AF2 


(A4) 


Analogous  to  the  continuous  case,  the  spectrum  is  obtained  by 
removing  the  delta  function  factor  enclosed  in  square  brackets. 


The  result  is  the  discrete  equivalent  of  (31),  showing  that  the 
above  prescriptions  generate  the  correct  random  forcing. 

The  simulations  of  Figs.  2  and  3  employed  the  following  pa¬ 
rameters:  N  =  256,  L  =  0.5  m,  and  At  =  43.2  s.  With  these 
choices,  the  highest  wave  number  ttN/L  is  about  1600  rad/m 
with  a  corresponding  decay  time  of  about  2000  s,  showing 
that  the  chosen  value  of  At  allows  more  than  adequate  sam¬ 
pling  of  the  decay  or  growth  of  the  smallest  features  included 
in  the  discrete  spectrum.  The  plots  of  Fig.  3  only  extend  to 
one-half  the  maximum  wave  number  in  order  to  clearly  display 
the  low-wave-number  parts  of  the  spectra.  Low-frequency  de¬ 
tails  are  somewhat  smeared  by  a  smoothing  operation  used  to 
reduce  statistical  noise  in  the  right-hand  side  plots  of  Fig.  3.  This 
filter  was  a  5  x  5  finite-impulse  response  filter  with  a  Gaussian 
shape  having  response  at  the  four  corners  smaller  than  the  center 
response  by  a  factor  of  6.  The  low  wave  numbers  of  interest 
are  of  the  order  of  the  cutoff  parameter  Kq  =  62.8  rad/ni.  The 
corresponding  decay  time  is  1.27  x  10^  s.  Thus,  the  simula¬ 
tion  must  run  for  more  than  3  x  10"^  steps  to  follow  the  growth 
of  the  largest  features.  Fortunately,  the  simulation  is  rather  ef¬ 
ficient  numerically.  Unlike  the  split-step  approximation  for  the 
parabolic  equation,  there  is  no  need  to  perform  Fourier  trans¬ 
forms  at  each  step.  The  discrete  Fourier  transform  Fi(7n,  71)  can 
be  marched  forward  in  time,  and  the  relief  function  in  coordi¬ 
nate  space  can  be  obtained  by  a  transform  occasionally  as  de¬ 
sired.  The  diffusion  propagator  matrix  cxp{-DK^^At)  need 
be  computed  only  once,  as  it  is  identical  in  each  time  step. 
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